Creating p-wave superfluids and topological excitations in optical lattices 
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We propose to realize a p-wave superfluid using bosons mixed with a single species of fermions in 
a deep optical lattice. We analyze with a self- consistent method its excitation spectrum in presence 
of a vortex, and we point out the range of interaction strengths in which the zero-energy mode with 
topological character exists on a finite optical lattice. Lattice effects are strongest close to fermionic 
half- filling: here the linearity of the low- lying spectrum is lost, and a new class of extended zero- 
energy modes with checkerboard structure and d-wave symmetry appears. 
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Pairing physics plays a central role in systems as di- 
verse as neutron stars, liquid Helium, solid state super- 
conductors and ultracold dilute gases [1]. Important 
advances in our understanding of this fascinating phe- 
nomenon came from the recent theoretical and exper- 
imental progresses made with two-component fermionic 
gases (2|,l3|]. In these systems, pairing occurs mainly in the 
s-wave channel and pairs have vanishing angular momen- 
tum. Ultracold atoms may provide the highly-controlled 
experimental framework needed to prepare and study 
unconventional superfluids, composed of pairs having a 
non-zero angular momentum. In two dimensions (2D), 
these superfluids feature anyonic excitations with a linear 
(Dirac-like) spectrum, and zero-energy excitations with 
topological character and non-Abelian statistics [4-7J. 
For p-wave symmetry, these excitations are analogous to 
the ones found in the v = b/2 Fractional Quantum Hall 
physics. Due to their insensibility to local perturbations, 
non-Abelian anyons are very appealing for applications 
in Topologically-protected Quantum Computation jsj. 

Although strong interactions between identical 
fermions have been obtained close to a p-wave Feshbach 
resonance (9|-11|, experiments so far failed to reach the 
superfluid regime since three-body collisions severely 
limit the lifetime of the gas. Recent theoretical pro- 
posals to realize long lived p-wave superfluids relied 
either on ensembles of 2D fermionic polar molecules 
dressed with radio-frequency |121] or on suppression of 
three-body collisions through the Quantum Zeno effect 
[l3]. Realization of p-wave and d-wave superfluids in 
shallow lattices has been proposed in Refs. (l3 . [lij . 
Large p-wave interactions are also predicted to arise 
between 3D lattice fermions in contact with 2D dipolar 
bosons [16]. 

In this paper, we propose to realize a stable p-wave 
superfluid in an optical lattice by doping a bosonic Mott 
state with a single fermionic species. We study its ex- 
citation spectrum in presence of a vortex by means of 
a self-consistent Bogolubov-de Gennes (BdG) approach. 
As shown in the seminal paper of Read and Green 




Figure 1: The dispersion E{\^)/t of a vortex- free chiral p- 
wave superfluid on a lattice, Eq. (|4]), vanishes linearly around 
the Dirac points. These gapless points appear, from left to 
right, at /i = — 4t, 0, 4t, and they separate different topological 
phases, each identified by a set {Ck}- The four indices of the 
set are or 1, depending on whether the momentum state 
with = {(0, 0), (tt, 0), (0, tt), (tt, tt)} is empty or occupied 
[25]. While the topological transition at /i = — 4t is analogous 
to the one found in the continuum, the ones at /i = 0, At are 
entirely due to lattice physics. 



a 2D p-wave superfluid in the continuum has a zero en- 
ergy excitation (Majorana mode) with exotic topological 
properties. In a lattice of finite extension, comparable to 
the ones found in current experiments, we find that the 
zero-energy Majorana mode survives only for intermedi- 
ate interaction strengths. Moreover, the lattice brings 
along a number of peculiar effects. In particular, we 
find that at half filling the spectrum becomes gapless and 
the gas undergoes a Topological Phase Transition (TFT). 
Close to criticality, a new class of extended zero-energy 
modes appears, characterized by a distinctive checker- 
board structure and d-wave symmetry. 

Our proposal to form a p-wave gas employs a Bose- 
Fermi mixture in a deep lattice. We choose the chem- 
ical potential for the bosons such that, in the absence 
of fermions, the bosons form a Mott state with unit fill- 
ing. As described in Ref. [17], adding identical fermions 
with filling < F < 1 yields a system that may 
be described by an effective Hamiltonian for composite 
fermions (fermion+boson, or fermion+bosonic hole, de- 
pending on the sign of the interspecies interaction). At 
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low temperatures the identical composite fermions can- 
not collide in s-wave due to Pauli exclusion principle, and 
the interaction has a dominant p-wave character. Su- 
perfluidity occurs when both Bose-Fermi and Bose-Bose 
interactions are repulsive and satisfy \Ubf\ > \Ubb\ > 0, 
giving an attractive (superexchange) nearest-neighbor in- 
teraction between fermion-hole composites; this interac- 
tion is characterized by U/t = 2{Ubf /Ubb — for the 
coupling constant between nearest neighbors of a Fermi- 
Hubbard model [l7|. Here t = 2J'^/Ubf is the tunneling 
energy of the composites, and J is the corresponding en- 
ergy for single particles. We assume for simplicity equal 
values of J for both bosons and fermions, but the scheme 
can be easily generalized to the case J5 ^ J/. The full 
tight-binding Hamiltonian for composites with chemical 
potential /i on a square lattice with N = L^L sites reads 

N u 

H = -fi^cla - t c\cj - ^ XI (1) 
i=i <i,j> <i,j> 

In order for this single-band approximation to be valid, 
the interband gap should be larger than both the ther- 
mal and interaction energies. This can easily achieved 
in typical deep lattice setups. Employing an s-wave Fes- 
hbach resonance to increase Ubf or decrease Ubb allows 
one to reach the large values of U/t needed to obtain, 
as we will see, bulk gaps of order the hopping t [18], 
and therefore (within a mean-field calculation) estimated 
critical temperatures of order the super-exchange inter- 
action [19||. Super-exchange effects have been recently 
observed in a bosonic gas [20], and with the develop- 
ment of novel cooling strategies (see e.g.0|) such physics 
seems to be at reach in the near future in lattice mix- 
tures. The desired control over the relevant parameters 
could for example be obtained in a ^^^Yb-^^^Yb mix- 
ture, where at zero magnetic field one has Ubf/Ubb ^ 4 
[22], or in a ^Li-^Li mixture, where U/t may be tuned 
employing well-characterized homonuclear and heteronu- 
clear s-wave Feshbach resonances [23]. It should be no- 
ticed that in our proposal three-body losses are expected 
to be severely limited, since we employ only s-wave Fesh- 
bach resonances, and in real (virtual) processes sites are 
never occupied by more than one (two) particle(s). 

Following the standard BCS treatment, we introduce 
the p-wave gap function Aij = U (ciCj), which is naturally 
antisymmetric due to the anticommuting properties of 
the fermionic operators q. Keeping only quadratic terms 
in the fluctuations around the mean-field, the Hamil- 
tonian may be diagonalized by means of the operators 
7n = Un{i)ci-\-Vn(i)cl {u = 1, . . . , 2A^), leading to the 
BdG equations (c.f. [6|): 

(-Alt') (:;)=-"(:■)■ 

Here hoUn{i) = -iiUnii) - '^n(j) is the ideal-gas 

tight-binding Hamiltonian, while the pairing operator A 



satisfies A''' = —A*, and acts on the quasiparticle wave- 
functions as Aun{i) = j) \jUn{j) (the two sums 
above run over j only). The 2N allowed energies appear 
in pairs of opposite sign since the Hamiltonian is symmet- 
ric under the exchange of {En^ipn} with {— <Ji?/^*)}, 
where ai is the first Pauli matrix and i/jn = {un^Vn)^ • 
Anyway, only the N excitations with En > are physical 
and will be considered in the following. As a consequence 
of this symmetry, if there exists a solution with ^0 = 0, 
the associated wavefunction may be chosen to satisfy the 
condition uq = Vq. The self-consistency condition, or 
gap-equation, reads 

A,j = V{c,c,) = UJ2 <(i)«n(i)tanh , (3) 

with kB the Boltzmann constant and T the temperature. 
The total number of particles is A/tot = — 
T.E^>^,i [M)? + f{E^){\un{i)? - M)?)] Where 
f{E) is the Fermi distribution, and the lattice filling is 
F = ATtot/AT. 

On a lattice, the dispersion relation for a homogeneous 
(vortex- free) chiral {px ± ipy) superfluid with order pa- 
rameter A(k) = Ah[s'm{kxd) ± ism{kyd)] reads (l^l 

^(k) = ^Je^\^h\^[^in\kxd)^sm\kyd)] (4) 

where ^ = —fi — 2t[cos{kxd) + cos{kyd)], and d is the lat- 
tice constant. This dispersion presents gapless points for 
three values of /i, i.e. ja = —4:t^ 0, 4t. In the absence of 
interactions, these values correspond respectively to av- 
erage fillings F = 0, 1/2, 1. When ^ the dispersion 
vanishes linearly around the gapless points, see Fig. [TJ 
i.e, it coincides with the one of the Dirac equation. 

Topological modes are known to arise as excitations 
bound to the core of a vortex in the order parameter. 
To model the vortex, we use open boundary conditions 
and solve self-consistently Eqs. (|2][3]) by taking the initial 
Ansatz for the gap Aij = xifi^'"^^' + /^•e^^^^]/2. Here 
fi = Ti/ \/rl ^ 2g^ is the amplitude of the gap around a 
vortex with core size a, and Vi and di are the polar co- 
ordinates of site i with respect to the core. The chirality 
factor X for Px^'^Vy pairing takes the values x = il when 
j = z ± X, and x = ii when j = z ± y. Finally, il' = ±1 
is the direction of rotation of the vortex with respect to 
the chirality. The results presented here are obtained on 
a lattice with 25 * 25 sites, which is experimentally feasi- 
ble, and large enough to obtain a satisfactory qualitative 
description of larger lattices [26]. 

In the absence of a confining potential, the physics in 
the lattice exhibits particle-hole symmetry around the 
half- filling point /i = 0. For /i < (/i > 0), one has 
< F < 1/2 (1/2 < F < 1), and the superfluid may be 
identified with a condensate of pairs of particles (holes). 
The gap function and the energy of the eigenstates de- 
pend only on the modulus of /i, but the wavefunctions u 
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Figure 2: Bulk gap Ao and lowest energy eigenvalues Ei for 
the case of a vortex rotating together with the chirality {w — 
+ 1). Energies in units of t. In Fig. [2]3, the thick line at 
l_x = -4t marks the TPT from strong 0000 to weak 1000 
pairing. The non-trivial TPT from weak 1000 to weak 1110 
pairing lies at F=0.5. In the lower figures (where U = 5t) 
the arrow marks the crossing of two eigenvalues causing the 
abrupt spreading of ipo close to the TCP at F = 1/2. The 
thin Unes in (d) represent the case w = —1. 



and V depend also on its sign, since the two phases 1000 
and 1110 are topologically distinguishable [25]. 

In the continuum, zero energy modes appear when the 
chemical potential becomes positive, the so-called weak 
pairing phase. Correspondingly, in a lattice the weak 
pairing phases are found for chemical potentials in the 
range — 4t < /i < and < /i < 4t for particles and holes, 
respectively. As shown in Fig. [21 at intermediate fillings 
and for sufficiently strong interactions (t < U < lOt) 
the T=0 bulk gap Aq (defined as the modulus of A^^ 
away from the vortex core) rises to values ~ t indicating 
large pairing correlations, and the low energy spectrum 
of the BdG Hamiltonian becomes linear, \En\ ^ ncuo, 
where n = 0, 1, 2, . . . (27|]. In particular, there exists a 
quasi-Majorana fermion excitation i/jq with uq Vq and 
an extremely small eigenvalue, Eq <^ Aq. For very strong 
attraction {U > lOt) the vortex core size shrinks below 
a lattice constant, and the lowest energy excitation loses 
its topological quasi-Majorana character. Moreover, too 
large values oiU/t should be avoided since they lead to 
phase separation [14, 28]. As intuitively expected from s- 
wave BCS theory, we have checked that a quasi-Majorana 
mode is present as long as the temperature is smaller than 
Ao, the bulk gap at T = 0. 

In Fig. [3] we show a typical ?/^o in the weak pairing 
phase; similarly to the continuum zero- modes [4-6], it ex- 
hibits oscillations which are exponentially damped away 
from the vortex core. When the vortex rotates in the 
(opposite to the) direction of the chirality, i.e., w = ^1 
{w = —1), IV^ol has a node (a maximum) at the vor- 
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Figure 3: Wavefunction uq — Vq at F ^ 1/4, w = —1, and 
U = 5t. Left: the modulus is exponentially localized around 
the vortex core. Right: the phase plot evidences the shape of 
the nodes adapting to the lattice structure. 
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tex core. The lattice effects are clearly reflected in four 
typical features: the energy of quasi-Majorana mode Eq 
exhibits oscillations as a function of F (see Fig.[2]l), cor- 
responding to the appearance of new nodes within the 
finite extension of the lattice; the nodal structure of the 
wavefunctions is aligned along the lattice axes to mini- 
mize energy; for the level spacing we find cjq ~ Aq, while 
in the continuum one has ujq ~ Aq/cf with cf the Fermi 
energy; finally, at half-filling the momentum distribution 
?/^o(k) is peaked at the minima (Dirac points) of the dis- 
persion relation of Fig. [TJcenter). 

The lattice effects are most strongly noticeable at half- 
filling (/i = 0), where the lattice spacing gets to coincide 
with half the period of the lowest energy excitation os- 
cillations. Already in the homogeneous, vortex free case, 
there exist two zero-energy "stripe modes", since the dis- 
persion vanishes at the gapless points k = (0, ±7r) and 
k = (±7r,0), see Eq. (j4]) and Fig. [TJcenter), and these 
modes may combine to form a checkerboard d-wave pat- 
tern. In presence of a vortex Eq. (j4]) does not hold, and 
the situation is more complex. The system minimizes its 
energy by arranging itself in a d-wave checkerboard pat- 
tern, the vortex core being occupied, see Fig. HI With 
open boundary conditions, however, this pattern is com- 
mensurate with the lattice only when L (the number of 
sites per lattice side) is odd: in this case, the correspond- 
ing mode has a vanishingly small eigenvalue [25]. If L is 
even, this zero is lifted and the spectrum contains two 
identical lowest eigenvalues. This situation is reminis- 
cent of the spectrum of a non-interacting chain, where 
the presence of a zero eigenvalue depends on its even or 
odd length. Lattice effects at /i = also destroy the lin- 
earity of the low energy part of the spectrum, see Fig. [2b. 

For odd L, ipo strongly depends on the direction of ro- 
tation of the vortex with respect to the chirality, as shown 
in Fig. m When the two counter-rotate {w = —1), the 
even symmetry of the wavefunction is compatible with 
the requirement of d-wave symmetry at half filling. As 
F approaches 1/2, i/jq smoothly evolves towards a wave- 
function tightly localized around the core. When instead 

= 1, as discussed above, away from half- filling the 
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Figure 4: Wavefunction i^o — '^^o ^it half-filling (/i = 0). Left: 
sketch of the d-wave checkerboard structure (+ and - refer to 
the sign of the Re(i/^o), while Im(i^o) = 0). Center and right: 
modulus of uo for w = —1 and w = respectively. 



lowest energy excitation has a node in the vortex core. 
This is not compatible with the checkerboard structure 
minimizing the energy at half-filling. As a consequence, 
approaching the Topological Critical Point (TCP) /i = 
the symmetry of the lowest energy excitation suddenly 
changes, and i/jq spreads over the whole lattice. The sud- 
den spreading of tjjo is signalled by the sharp crossing 
of the two lowest eigenvalues, as marked by an arrow 
in Fig. [2]3-d. On a finite lattice, the detection of this 
extended state would therefore be a clear experimental 
signature of the vicinity of the TCP. 

The general features discussed above persist even in 
the presence of an external harmonic confinement, pro- 
vided this varies sufficiently slowly in the central part of 
the lattice. In particular, we find that the low-energy 
section of the (positive) spectrum in the weak pairing 
phase is still approximately linear and admits an eigen- 
value Eq <C Aq. Nonetheless, the confinement destroys 
the aforementioned particle-hole symmetry, and the spa- 
tial inhomogeneity makes it difficult to observe the for- 
mation of the checkerboard state around half-filling. To 
avoid this problem, optical potentials formed by higher 
order Laguer re- Gauss beams, or by "light sheets", could 
be used; they yield potentials of the form oc (x^^ + y'^^)^ 
characterized by a rather uniform central region and 
steep confining walls. Simulations for the experimen- 
tally feasible case u = 2, w = 1, and F 0.22, lead to 
/i = 0, and exhibit quite nicely a checkerboard quasi-zero 
mode with significant extension. This structure could 
be detected experimentally via spatially-resolved polar- 
ization spectroscopy [29], while Majorana modes may be 
observed via radio- frequency absorption jsoj- 

In conclusion, we have proposed a novel scheme to re- 
alize stable p-wave superfluids in optical lattices, and we 
have shown that even in presence of a relatively small 
lattice, the lowest energy excitation preserves its topo- 
logical character. Optimal filling factors and interaction 
strengths to obtain Majorana modes are F ^ 1/4 (or 
equivalently F ~ 3/4) and U/t ~ 5, where the ratio 
Eq/Ao is smallest. At the critical point fJ^ = 0, the sys- 
tem still possesses a zero mode, but the linear structure 
of the low energy part of the spectrum is lost. Here, the 
system favors the formation of checkerboard order. If the 



vortex rotates in the same direction as the chirality, the 
vicinity of the topological phase transition is signalled 
by the abrupt spreading of the lowest energy excitation 
wavefunction. These results are robust for temperatures 
below the super-exchange interaction and even in pres- 
ence of an external confinement. 
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